#########Estancona and Tiscornia 
#########Replication for `From Cocaine to Avocados`

rm(list=ls(all=TRUE))
library(tidyverse)
library(foreign)
library(lme4)
library(readstata13)
library(DataCombine)
library(ggplot2)
library(stargazer)
library(dplyr)
library(tidyr)
library(haven)
library(foreign)
library(reshape2)
library(tibbletime)
library(stringi)
library(fixest)
library(modelsummary)
library(interflex)
library(vdemdata)
library(miceadds)
library(ggrepel)
library(caTools)
library(gridExtra)
library(ggpubr)
library(xtable)
library(scales)


setwd("")

####################bring in data
####Mexico files
agri<-read.csv("Agricultural_Production.csv")
mexdat<-read.csv("mexdat_IO.csv")
toast<-read.csv("Avo_Toast_Searches.csv")
toast_us<-read.csv("Avo_Toast_Searches_US.csv")
corn<-read.csv("corn.csv")
strawberries<-read.csv("strawberries.csv")
limes<-read.csv("limes.csv")

####general data files
homicide_total<-read.csv("homicide_total_rate_and_count_2019.csv")
xc_data<-read.csv("cross_national_data.csv")

#############
#######################################
#plots and tables to produce, in order in the paper 
#######################################
#############

#################################
#################################For Mexico section:

######mexico case plot (situating case) - mexico_global_scatterplot_2.pdf
#Figure 1 

g_av<-xc_data[which(xc_data$year>=1993),]
g_av<-g_av[,-1]
g_av<-g_av[, c(1,2,9,22)]

homicide_rates<-homicide_total[which(homicide_total$Indicator=="Homicide rate"),]

g_av<-merge(g_av, homicide_rates, by=c("Country", "year"))
g_av<-g_av[,c(1,2,3,4,7)]
colnames(g_av)[5]<-"rate"

g_av<-unique(g_av)

g_av<- g_av %>% 
  group_by(Country) %>%
  mutate(av_ag=mean(export_av, na.rm=T),hom_av=mean(rate, na.rm=T))

g_av_2<-g_av[which(g_av$hom_av>=10),]

g_av<-g_av[, c(1,6,7)]

g_av<-unique(g_av)

mexico_scatterplot_all<-ggplot(g_av, aes(x=av_ag, y=hom_av)) + 
  geom_point()+
  geom_point(data=g_av[g_av$Country == "Mexico",],color="red",size=3)+
  geom_label_repel( 
    data=g_av %>% filter(hom_av>=10 & Country!="Mexico"),
    aes(label=Country))+
  geom_label( 
    data=g_av %>% filter(Country=="Mexico"),
    aes(label=Country),
    label.size=.1,
    vjust=1.2,
    color="red")+
  xlab("Average Agricultural Export Value")+
  ylab("Average Homicides")+
  theme(aspect.ratio=0.75)+
  theme_bw()

mexico_scatterplot_all

pdf("mexico_global_scatterplot_2.pdf")
mexico_scatterplot_all
dev.off()



######avocado production and presence plot - images/crim-orgs-map-gray.png
#Figure 2
#see additional QGIS files for production

#Change in Avocado Export Value Share, Criminal Threat, and Homicides - Avos_tab
#Table 3

test_years_avo_RR_3<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=mexdat)
summary(test_years_avo_RR_3)

test_years_avo_RR_no_cont<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs, cluster=~cve_inegi, data=mexdat)
summary(test_years_avo_RR_no_cont)

#reporting
modelsummary(list(test_years_avo_RR_no_cont, test_years_avo_RR_3), fmt=3, stars=T, output="latex")


####marginal_effect_Mex_no_bin.pdf - marginal effects
#Figure 3

t4_no_binning_no_fe <- interflex(Y = "homs_led_no_NA_std", D = "ev_share_change_std", X = "Number_Orgs", Z = c("election_current_year", "finpub_ing_brutos_std", "pob_total_est_std", "seg_agentes_fc_std"), 
                                 D.ref=max(na.omit(mexdat$ev_share_change_std)),
                                 data = as.data.frame(mexdat), estimator = "linear", vcov.type="cluster",
                                 cl="cve_inegi", na.rm=TRUE, ylab="Marginal Effect",
                                 xlab="Criminal Threat", xlim=c(0, 8.75), theme.bw=T)
                                 file="/marginal_effect_Mex_no_bin.pdf")

t4_no_binning_no_fe


#prediction_plot_Mex_norgs.pdf - predicted probs
#Figure 4

test_years_avo_RR_4<-feols(homs_led_no_NA~ev_share_change+Number_Orgs+ev_share_change*Number_Orgs+election_current_year+finpub_ing_brutos+pob_total_est+seg_agentes_fc, cluster=~cve_inegi, data=mexdat)
summary(test_years_avo_RR_4)


x_vec_ev_change<-seq(0, max(na.omit(mexdat$ev_share_change)), length.out=5000)


new_dat_norgs_1<-data.frame(x_vec_ev_change, 0,
                            0, median(na.omit(mexdat$finpub_ing_brutos)), median(na.omit(mexdat$pob_total_est)), 
                            median(na.omit(mexdat$seg_agentes_fc)))

new_dat_norgs_2<-data.frame(x_vec_ev_change, 1,
                            0, median(na.omit(mexdat$finpub_ing_brutos)), median(na.omit(mexdat$pob_total_est)), 
                            median(na.omit(mexdat$seg_agentes_fc)))


new_dat_norgs_3<-data.frame(x_vec_ev_change,9,
                            0, median(na.omit(mexdat$finpub_ing_brutos)), median(na.omit(mexdat$pob_total_est)), 
                            median(na.omit(mexdat$seg_agentes_fc)))

colnames(new_dat_norgs_1)<-c("ev_share_change", "Number_Orgs", "election_current_year", "finpub_ing_brutos", "pob_total_est", "seg_agentes_fc") 
colnames(new_dat_norgs_2)<-c("ev_share_change", "Number_Orgs", "election_current_year", "finpub_ing_brutos", "pob_total_est", "seg_agentes_fc") 
colnames(new_dat_norgs_3)<-c("ev_share_change", "Number_Orgs", "election_current_year", "finpub_ing_brutos", "pob_total_est", "seg_agentes_fc") 

fit_norgs_1.2<-predict(test_years_avo_RR_4, newdata=new_dat_norgs_1, interval="confidence")

fit_norgs_1.3<-predict(test_years_avo_RR_4, newdata=new_dat_norgs_2, interval="confidence")

fit_norgs_1.4<-predict(test_years_avo_RR_4, newdata=new_dat_norgs_3, interval="confidence")

plot_dat_norgs_1.2<-data.frame(fit_norgs_1.2, fit_norgs_1.3, fit_norgs_1.4, new_dat_norgs_1)


pdf("/Users/chelseaestancona/Dropbox/Chelsea and Lucia Projects/Criminal Groups Paper/Data/prediction_plot_Mex_norgs.pdf")
ggplot(data=plot_dat_norgs_1.2)+ 
  labs(x="Change in Share of Avocado Export Value", y="Homicide Count")+
  geom_line(aes(x=ev_share_change, y=fit_norgs_1.2[,1], color="No Criminal Threat"),linetype="solid")+
  geom_line(aes(x=ev_share_change, y=fit_norgs_1.4[,1], color="Max Criminal Threat"), linetype="solid")+
  geom_ribbon(aes(x=ev_share_change, ymin=fit_norgs_1.2[,3], ymax=fit_norgs_1.2[,4]), alpha=.3, fill="royalblue")+
  geom_ribbon(aes(x=ev_share_change, ymin=fit_norgs_1.4[,3], ymax=fit_norgs_1.4[,4]), alpha=.3, fill="firebrick1")+
  theme_set(theme_bw())+
  theme(legend.title=element_text(size=20), legend.position="bottom", panel.grid.minor.x=element_blank(),  panel.grid.minor.y=element_blank())+
  scale_fill_manual(values=c("royalblue", "firebrick"), name="")+
  scale_color_manual(values=c("firebrick", "royalblue"), name="")+
  scale_x_continuous(labels = label_percent())
dev.off()


####parallel trends -pt_plot.pdf
#Figure 5

mdf_avos<-mexdat[which(mexdat$avos_2==1),]

mdf_no_avos<-mexdat[which(mexdat$avos_2==0),]

avos_no_crim<-mdf_avos[which(mdf_avos$Number_Orgs==0),]
avos_crim<-mdf_avos[which(mdf_avos$Number_Orgs>=1),]
no_avos_no_crim<-mdf_no_avos[which(mdf_no_avos$Number_Orgs==0),]
no_avos_crim<-mdf_no_avos[which(mdf_no_avos$Number_Orgs>=1),]

#avocados no crime
anc_h<-avos_no_crim[,c("cve_inegi", "year", "homs_led_no_NA", "ev_share_change")]
anc_h<-anc_h %>%
  group_by(year) %>%
  mutate(av_homs=mean(na.omit(homs_led_no_NA)), av_ev_share=mean(na.omit(ev_share_change)))

anc_h<-anc_h[,c("year", "av_homs", "av_ev_share")]
anc_h<-unique(anc_h)

#avocados, crime
ac_h<-avos_crim[,c("cve_inegi", "year", "homs_led_no_NA", "ev_share_change")]
ac_h<-ac_h %>%
  group_by(year) %>%
  mutate(av_homs=mean(na.omit(homs_led_no_NA)), av_ev_share=mean(na.omit(ev_share_change)))

ac_h<-ac_h[,c("year", "av_homs", "av_ev_share")]
ac_h<-unique(ac_h)

#no avocados, no crime
nanc_h<-no_avos_no_crim[,c("cve_inegi", "year", "homs_led_no_NA", "ev_share_change")]
nanc_h<-nanc_h %>%
  group_by(year) %>%
  mutate(av_homs=mean(na.omit(homs_led_no_NA)), av_ev_share=mean(na.omit(ev_share_change)))

nanc_h<-nanc_h[,c("year", "av_homs", "av_ev_share")]
nanc_h<-unique(nanc_h)

#no avocados, crime
nac_h<-no_avos_crim[,c("cve_inegi", "year", "homs_led_no_NA", "ev_share_change")]
nac_h<-nac_h %>%
  group_by(year) %>%
  mutate(av_homs=mean(na.omit(homs_led_no_NA)), av_ev_share=mean(na.omit(ev_share_change)))

nac_h<-nac_h[,c("year", "av_homs", "av_ev_share")]
nac_h<-unique(nac_h)

colnames(anc_h)<-c("year", "anc_h", "anc_ev_share")
colnames(ac_h)<-c("year", "ac_h", "ac_ev_share")
colnames(nanc_h)<-c("year", "nanc_h", "nanc_ev_share")
colnames(nac_h)<-c("year", "nac_h", "nac_ev_share")

pt<-merge(anc_h, ac_h, by="year")
pt<-merge(pt, nanc_h, by="year")
pt<-merge(pt, nac_h, by="year")


pdf("/pt_plot.pdf")
ggplot(pt, aes(x=year, y=anc_h, color="Avocados, No Criminal Presence")) + geom_line(group=1) + geom_point()+
  geom_point(aes(y=ac_h, color="Avocados, Criminal Presence"))+ geom_line(aes(y=ac_h, group=1, color="Avocados, Criminal Presence"))+
  geom_point(aes(y=nac_h, color="No Avocados, Criminal Presence"))+ geom_line(aes(y=nac_h, group=1, color="No Avocados, Criminal Presence"))+
  geom_point(aes(y=nanc_h, color="No Avocados, No Criminal Presence"))+ geom_line(aes(y=nanc_h, group=1, color="No Avocados, No Criminal Presence"))+
  geom_vline(xintercept=2006, linetype="dotted", color="red")+
  geom_text(aes(x=2006, label="\nAvocado Toast", y=40), colour="red", angle=90)+
  labs(x="Year", y="Average Municipal Homicides", color="Type of Municipality")
dev.off()

####Homicides, Avocado Toast Search Popularity and Criminal Threat - toast_tab
#Table 4


mdf_avos$increase_search<-ifelse(mdf_avos$change_search>0,1,0)
mdf_no_avos$increase_search<-ifelse(mdf_no_avos$change_search>0,1,0)

test_avo_toast_controls<-feols(homs_led_no_NA_std~prop_years_2004+increase_search+increase_search*prop_years_2004+
                                 election_current_year+finpub_ing_brutos_std+pob_total_est_std+
                                 seg_agentes_fc_std+as.factor(year),
                               cluster=~cve_inegi,
                               data=mdf_avos)

summary(test_avo_toast_controls)


test_no_avo_toast_controls<-feols(homs_led_no_NA_std~prop_years_2004+increase_search+increase_search*prop_years_2004+
                                    election_current_year+finpub_ing_brutos_std+pob_total_est_std+
                                    seg_agentes_fc_std+as.factor(year),
                                  cluster=~cve_inegi,
                                  data=mdf_no_avos)

summary(test_no_avo_toast_controls)

modelsummary(list(test_avo_toast_controls, test_no_avo_toast_controls), fmt=3, stars=T, output="latex")


###avocado exports, toast - marginal_effect_avos.pdf
#Figure 6 a

avos_no_binning <- interflex(Y = "homs_led_no_NA_std", D = "increase_search", X = "prop_years_2004", Z = c("election_current_year", "finpub_ing_brutos_std", "pob_total_est_std", "seg_agentes_fc"), 
                             data = as.data.frame(mdf_avos), estimator = "linear", vcov.type="cluster",
                             cl="cve_inegi", na.rm=TRUE, ylab="Marginal Effect",
                             treat.type="discrete", 
                             xlab="Criminal Threat", ylim=c(0,1.5), xlim=c(0,.8))
                             file="/marginal_effect_avos.pdf")

avos_no_binning

###no avocado exports, toast - marginal_effect_no_avos.pdf
#Figure 6 b

no_avos_no_binning <- interflex(Y = "homs_led_no_NA_std", D = "increase_search", X = "prop_years_2004", Z = c("election_current_year", "finpub_ing_brutos_std", "pob_total_est_std", "seg_agentes_fc"), 
                                D.ref=c(0,1), base=0,   
                                data = as.data.frame(mdf_no_avos), estimator = "linear", vcov.type="cluster",
                                cl="cve_inegi", na.rm=TRUE, ylab="Marginal Effect",
                                treat.type="discrete", 
                                xlab="Criminal Threat", show.all=T, ylim=c(0,1.5), xlim=c(0,.8),
                                file="/marginal_effect_no_avos.pdf")

no_avos_no_binning

###avocado toast, predicted homicides - prediction_plot_avos_2_1.pdf
#Figure 7

x_vec_prop<-seq(min(na.omit(mdf_avos$prop_years_2004)), max(na.omit(mdf_avos$prop_years_2004)), length.out=5000)

test_avo_toast_controls_2_std<-feols(homs_led_no_NA~increase_search+prop_years_2004+increase_search*prop_years_2004+
                                       election_current_year+finpub_ing_brutos+pob_total_est+seg_agentes_fc+as.factor(year),
                                     cluster=~cve_inegi,
                                     data=mdf_avos)
summary(test_avo_toast_controls_2_std)

test_no_avo_toast_controls_2_std<-feols(homs_led_no_NA~increase_search+prop_years_2004+increase_search*prop_years_2004+
                                          election_current_year+finpub_ing_brutos+pob_total_est+seg_agentes_fc+as.factor(year),
                                        cluster=~cve_inegi,
                                        data=mdf_no_avos)
summary(test_no_avo_toast_controls_2_std)


preds2_avos<-plot_predictions(test_avo_toast_controls_2_std, 
                              condition=list("increase_search", "prop_years_2004"=c(0,0.0714285714285714, 0.785714285714286)), 
                              data=mdf_avos)+labs(x="Avocados", y="",
                                                  fill="Criminal Threat", colour="Criminal Threat")+
  scale_fill_manual(values=c('royalblue', 'springgreen4', 'firebrick'),labels=c("None", "Low", "High"))+
  scale_colour_manual(values=c('royalblue', 'springgreen4', 'firebrick'), labels=c("None", "Low", "High"))+ylim(-50,300)+theme_bw()

preds2_avos

preds2_no_avos<-plot_predictions(test_no_avo_toast_controls_2_std, 
                                 condition=list("increase_search", "prop_years_2004"=c(0,0.0714285714285714, 0.785714285714286)), 
                                 data=mdf_no_avos)+labs(x="No Avocados", y="", fill="Criminal Threat", colour="Criminal Threat")+
  scale_fill_manual(values=c('royalblue', 'springgreen4', 'firebrick'),labels=c("None", "Low", "High"))+
  scale_colour_manual(values=c('royalblue', 'springgreen4', 'firebrick'), labels=c("None", "Low", "High"))+ylim(-50,300)+theme_bw()

preds2_no_avos

c_avos_preds<-ggarrange(preds2_avos, preds2_no_avos, common.legend=T)

c_avos_preds_final<-annotate_figure(c_avos_preds,
                                    #top = text_grob("Visualizing len", color = "red", face = "bold", size = 14),
                                    bottom = text_grob("Increase in Avocado Toast Search Popularity", color = "black"),
                                    left = text_grob("Predicted Homicides", color = "black", rot = 90),
                                    #right = "I'm done, thanks :-)!",
                                    #fig.lab = "Figure 1", fig.lab.face = "bold"
)

pdf("/prediction_plot_avos_2_1.pdf")
c_avos_preds_final
dev.off()


####Corn, Strawberry, and Lime Production, Criminal Threat, and Homicides in non-Avocado Growing Municipalities- crops tab
#Table 5

#corn
test_years_corn_only_no_fe_2<-feols(homs_led_no_NA_std~pv_change_corn_std+Number_Orgs+pv_change_corn_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=corn_no_avos)
summary(test_years_corn_only_no_fe_2)

#strawberries
test_years_strawb_only_no_fe_2<-feols(homs_led_no_NA_std~pv_change_strawb_std+Number_Orgs+pv_change_strawb_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=strawb_no_avos)
summary(test_years_strawb_only_no_fe_2)

#limes
test_years_limes_only_no_fe_2<-feols(homs_led_no_NA_std~pv_change_limes_std+Number_Orgs+pv_change_limes_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=limes_no_avos)
summary(test_years_limes_only_no_fe_2)

#reporting
modelsummary(list(test_years_corn_only_no_fe_2, test_years_strawb_only_no_fe_2, test_years_limes_only_no_fe_2), fmt=3, stars=T, output="latex")


#######################################################
########################################General tests:

#####Change in Export Value Share, Criminal Threat, and Homicides, tab_main
## Table 6

RR_main_model_no_fe<-lm(led_homs_std~product_change_export_value_share+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                          sitc_eci_2+v2csanmvch_12+v2csanmvch_12*product_change_export_value_share+v2x_rule+v2xnp_regcorr, cluster=~Country, data=xc_data)

summary(RR_main_model_no_fe)

RR_main_model_no_fe_base<-lm(led_homs_std~product_change_export_value_share+
                               v2csanmvch_12+v2csanmvch_12*product_change_export_value_share, cluster=~Country, data=xc_data)

summary(RR_main_model_no_fe_base)

#reporting
modelsummary(list(RR_main_model_no_fe_base, RR_main_model_no_fe), fmt=3, stars=T, output="latex")

#######prediction_xc_RR.pdf, general predicted probabilities
##Figure 8
criminal_presence<-c(min(na.omit(xc_data$v2csanmvch_12)), mean(na.omit(xc_data$v2csanmvch_12)), quantile(na.omit(xc_data$v2csanmvch_12))[4], max(na.omit(xc_data$v2csanmvch_12)))
x_vec_a<-seq(min(na.omit(xc_data$product_change_export_value_share)), max(na.omit(xc_data$product_change_export_value_share)), length.out=5000)


#predictions for different values of criminal threat across value of export change
new_dat_box_1<-data.frame(criminal_presence[1], median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                          median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_a, 
                          median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))

new_dat_box_2<-data.frame(criminal_presence[2], median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                          median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_a, 
                          median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))

new_dat_box_3<-data.frame(criminal_presence[3], median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                          median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_a, 
                          median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))       

new_dat_box_4<-data.frame(criminal_presence[4], median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                          median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_a, 
                          median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))   


colnames(new_dat_box_1)<-c("v2csanmvch_12", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", 
                           "product_change_export_value_share", "v2x_rule", "v2xnp_regcorr")

colnames(new_dat_box_2)<-c("v2csanmvch_12", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", 
                           "product_change_export_value_share", "v2x_rule", "v2xnp_regcorr")

colnames(new_dat_box_3)<-c("v2csanmvch_12", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", 
                           "product_change_export_value_share", "v2x_rule", "v2xnp_regcorr")

colnames(new_dat_box_4)<-c("v2csanmvch_12", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", 
                           "product_change_export_value_share", "v2x_rule", "v2xnp_regcorr")

fit_r_1<-predict(RR_main_model_no_fe_no_std, newdata=new_dat_box_1, interval="confidence")


fit_r_2<-predict(RR_main_model_no_fe_no_std, newdata=new_dat_box_2, interval="confidence")


fit_r_3<-predict(RR_main_model_no_fe_no_std, newdata=new_dat_box_3, interval="confidence")


fit_r_4<-predict(RR_main_model_no_fe_no_std, newdata=new_dat_box_4, interval="confidence")


plot_dat_r<-data.frame(fit_r_1, fit_r_2, fit_r_3, fit_r_4, new_dat_box_1)
plot_dat_r<-plot_dat_r[which(plot_dat_r$fit>0 & plot_dat_r$fit.1>0 & plot_dat_r$fit.3>0),]

#pdf("prediction_xc_RR.pdf")
ggplot(plot_dat_r)+ 
  labs(x="Export Value Share Change", y="Predicted Homicides")+
  geom_line(aes(x=product_change_export_value_share, y=plot_dat_r[,1], color="No Criminal Threat"),linetype="solid")+
  geom_line(aes(x=product_change_export_value_share, y=plot_dat_r[,4], color="Mean Criminal Threat"),linetype="solid")+
  geom_line(aes(x=product_change_export_value_share, y=plot_dat_r[,10], color="Max Criminal Threat"),linetype="solid")+
  geom_ribbon(aes(x=product_change_export_value_share, ymin=plot_dat_r[,2], ymax=plot_dat_r[,3]), alpha=.3, fill="royalblue")+
  geom_ribbon(aes(x=product_change_export_value_share, ymin=plot_dat_r[,5], ymax=plot_dat_r[,6]), alpha=.3, fill="springgreen1")+
  geom_ribbon(aes(x=product_change_export_value_share, ymin=plot_dat_r[,11], ymax=plot_dat_r[,12]), alpha=.3, fill="firebrick1")+
  theme_set(theme_bw())+
  theme(legend.title=element_text(size=20), legend.position="bottom", panel.grid.minor.x=element_blank(),  panel.grid.minor.y=element_blank())+
  scale_fill_manual(values=c("firebrick","springgreen", "royalblue"), name="")+
  scale_color_manual(values=c("firebrick", "springgreen", "royalblue"), name="")
#dev.off()

######################################
#################################Appendix:

####Change in Avocado Export Value Share, Criminal Threat, and Homicides, Fixed Effects Included, Avocados_fe_tab
#Table 7

test_years_avo_RR<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs| year + cve_inegi, cluster=~cve_inegi, data=mexdat)
summary(test_years_avo_RR)


test_years_avo_RR_2<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std
                           | year + cve_inegi, cluster=~cve_inegi, data=mexdat)


#reporting
modelsummary(list(test_years_avo_RR, test_years_avo_RR_2), fmt=3, stars=T, output="latex")


#####Change in Avocado Export Value, Avocados_ev_tab
#Table 8

test_ev_change_no_fe<-feols(homs_led_no_NA_std~export_value_change+Number_Orgs+export_value_change*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=mexdat)
summary(test_ev_change_no_fe)

test_ev_change_no_fe_base<-feols(homs_led_no_NA_std~export_value_change+Number_Orgs+export_value_change*Number_Orgs, cluster=~cve_inegi, data=mexdat)
summary(test_ev_change_no_fe_base)

#reporting
modelsummary(list(test_ev_change_no_fe_base, test_ev_change_no_fe), fmt=3, stars=T, output="latex")

#####Change in Local Avocado Price, Avocados_lp_tab
#Table 9

test_lp_change_no_fe<-feols(homs_led_no_NA_std~local_price_change+Number_Orgs+local_price_change*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=mexdat)
summary(test_lp_change_no_fe) 

test_lp_change_no_fe_base<-feols(homs_led_no_NA_std~local_price_change+Number_Orgs+local_price_change*Number_Orgs, cluster=~cve_inegi, data=mexdat)
summary(test_lp_change_no_fe_base) 

modelsummary(list(test_lp_change_no_fe_base, test_lp_change_no_fe), fmt=3, stars=T, output="latex")

####Change in Avocado Export Value Share, Criminal Threat, and Homicides Scaled by Population, Avocados_rate_tab
#Table 10

mexdat$homs_led_no_NA_rate<-mexdat$homs_led_no_NA/mexdat$pob_total_est

ev_share_rate_3<-feols(homs_led_no_NA_rate~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+seg_agentes_fc_std, cluster=~cve_inegi, data=mexdat)
summary(ev_share_rate_3)

ev_share_rate_3_base<-feols(homs_led_no_NA_rate~ev_share_change_std+Number_Orgs+ev_share_change_std*Number_Orgs, cluster=~cve_inegi, data=mexdat)
summary(ev_share_rate_3_base)

#reporting
modelsummary(list(ev_share_rate_3, ev_share_rate_3_base), fmt=6, stars=T, output="latex")


###Change in Avocado Export Share and Criminal Threat, No Interaction Term
#Table 11

ev_share_no_int_3<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs+election_current_year+finpub_ing_brutos_std+seg_agentes_fc_std, cluster=~cve_inegi, data=mexdat)
summary(ev_share_no_int_3)

ev_share_no_int_3_base<-feols(homs_led_no_NA_std~ev_share_change_std+Number_Orgs, cluster=~cve_inegi, data=mexdat)
summary(ev_share_no_int_3_base)

#reporting
modelsummary(list(ev_share_no_int_3_base, ev_share_no_int_3), fmt=3, stars=T, output="latex")

###Homicides and Change in Corn Production Value During Price Jump, corn_price_shock_tab
#Table 12
corn_lim<-corn[which(corn$year>2004 & corn$year<2009),]

test_years_corn_only_lim_no_fe_2<-feols(homs_led_no_NA_std~pv_change_corn_std+Number_Orgs+pv_change_corn_std*Number_Orgs+election_current_year+finpub_ing_brutos_std+pob_total_est_std+seg_agentes_fc_std, cluster=~cve_inegi, data=corn_lim)
summary(test_years_corn_only_lim_no_fe_2)

modelsummary(test_years_corn_only_lim_no_fe_2, fmt=3, stars=T, output="latex")


###marginal_effect_Mex.pdf - binned ME
#Figure 9

RR_mex_me <- interflex(Y = "homs_led_no_NA", D = "ev_share_change", X = "Number_Orgs", Z = c("election_current_year", "finpub_ing_brutos", "pob_total_est", "seg_agentes_fc"), 
                       D.ref=max(na.omit(mexdat$ev_share_change)),
                       data = as.data.frame(mexdat), estimator = "binning", vcov.type="cluster",
                       cl="cve_inegi", na.rm=TRUE, ylab="Marginal Effect", theme.bw=T,
                       xlab="Criminal Threat", file="/marginal_effect_mex.pdf")

RR_mex_me

####density_avos_export_share_change.pdf, Density plot, export val share
#Figure 10

pdf("/density_avos_export_share_change.pdf")
ggplot(data=mexdat)+ 
  aes(x=ev_share_change)+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.9)+
  labs(x="Change in Export Value Share", y="Density")+
  scale_x_continuous(limits = c(-0.0001, 0.0001))+
  theme_bw()
dev.off()


####density_Mex_criminal_threat.pdf, Density plot, # of crim orgs
#Figure 11

pdf("/density_Mex_criminal_threat.pdf")
ggplot(data=mexdat)+ 
  aes(x=Number_Orgs)+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.9)+
  labs(x="Criminal Threat", y="Density")+
  theme_bw()
dev.off()


###Change in Avocado Export Value Share, Criminal Presence, and Homicides - appendix_avocados
#Table 13

test_prop_years_base<-feols(homs_led_no_NA_std~prop_years_2004+ev_share_change_std+ev_share_change_std*prop_years_2004+as.factor(year),
                            data=mexdat)
test_prop_years_base

test_prop_years<-feols(homs_led_no_NA_std~prop_years_2004+ev_share_change_std+ev_share_change_std*prop_years_2004+
                                    election_current_year+finpub_ing_brutos_std+pob_total_est_std+
                                    seg_agentes_fc_std+as.factor(year),
                                  data=mexdat)
test_prop_years


test_prop_years_no_fe<-feols(homs_led_no_NA_std~prop_years_2004+ev_share_change_std+ev_share_change_std*prop_years_2004+
                         election_current_year+finpub_ing_brutos_std+pob_total_est_std+
                         seg_agentes_fc_std,
                       data=mexdat)
test_prop_years_no_fe

modelsummary(list(test_prop_years_base, test_prop_years, test_prop_years_no_fe), omit=c("year", "Country"), fmt=3, stars=T, output="latex")


###`Avocado Toast' Search Popularity, Number of OCGs, and Homicides
#Table 14

test_avo_toast_norgs<-feols(homs_led_no_NA_std~Number_Orgs+change_search_std+change_search_std*Number_Orgs,
                               data=mdf_avos)
test_avo_toast_norgs

test_no_avo_toast_norgs<-feols(homs_led_no_NA_std~Number_Orgs+change_search_std+change_search_std*Number_Orgs,
                            data=mdf_no_avos)
test_no_avo_toast_norgs

modelsummary(list(test_avo_toast_norgs, test_no_avo_toast_norgs),fmt=3, stars=T, output="latex")

########Descriptive Statistics, Dep. Var. and Continuous Indep. Vars. (Global)
########Descriptive Statistics, Discrete Independent Variables, Alternative Tests (Global)

#Table 15

test_for_vars<-lm(led_homs_std~shock_ev_share_any+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                  sitc_eci_2+v2csanmvch_12+v2csanmvch_12*shock_ev_share_any+v2x_rule+v2xnp_regcorr, data=xc_data)

test_for_vars_no_standardized<-lm(homs_led~shock_ev_share_any+Population+GDPPC+Ag_Emp_Percent+conflict_dummy+
                                  sitc_eci_2+v2csanmvch_12+v2csanmvch_12*shock_ev_share_any+v2x_rule+v2xnp_regcorr, data=xc_data)


summarydat<-model.frame(test_for_vars)
continuous<-as.data.frame(rbind(summary(summarydat$pop_std), summary(summarydat$gdppc_std), summary(summarydat$Ag_Emp_Percent_std),
                                summary(summarydat$sitc_eci_2), summary(summarydat$v2csanmvch_12), summary(summarydat$v2x_rule), summary(summarydat$v2xnp_regcorr), summary(summarydat$led_homs_std)))
rownames(continuous)<-c("Population (Standardized)", 
                        "GDPPC (Standardized)", "Percent Employed in Agriculture", "ECI", "Criminal Presence (Av. Coder)",
                        "Rule of Law", "Corruption", "Homicide Count (Standardized)")
xtable(continuous[,c(1,3,4,6)]) #limit to needed columns 

#Table 18
discrete<-as.data.frame(rbind(table(summarydat$conflict_dummy), table(summarydat$shock_ev_share_any)))
rownames(discrete)<-c("Conflict", "Shock")
colnames(discrete)<-c("No", "Yes")
xtable(discrete)

#######marginal_effect_xc.pdf- global ME 
#Figure 12 

RR_mod_me <- interflex(Y = "homs_led", D = "product_change_export_value_share", X = "v2csanmvch_12", Z = c("Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", "v2x_rule", "v2xnp_regcorr"), 
                       D.ref=max(na.omit(xc_data$product_change_export_value_share)),
                       data = as.data.frame(xc_data), estimator = "binning", vcov.type="cluster",
                       cl="Country", na.rm=TRUE, ylab="Marginal Effect", cutoffs=c(0.25,0.75), theme.bw=T,
                       xlab="Criminal Threat", file="")

RR_mod_me

######density_export_share.pdf - dist. of change in EV share
#Figure 13
pdf("/density_export_share_change.pdf")
ggplot(data=xc_data)+ 
  aes(x=product_change_export_value_share)+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.9)+
  labs(x="Product Level Change in Export Value Share", y="Density")+
  scale_x_continuous(limits = c(-0.05, 0.05))+
  scale_y_continuous(limits = c(0, 6250))+
  theme_bw()
dev.off()

######density_criminal_threat.pdf - dist. of criminal threat (global)
#Figure 14
pdf("/density_criminal_threat.pdf")
ggplot(data=xc_data)+ 
  aes(x=v2csanmvch_12)+
  geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.9)+
  labs(x="Criminal Threat", y="Density")+
  scale_x_continuous(limits = c(0,max(na.omit(g18$v2csanmvch_12))))+
  #scale_y_continuous(limits = c(0, 6250))
  theme_bw()
dev.off()

#######Change in Export Value Share, Criminal Threat, and Homicides, no Interaction Term - xc_change_no_int
##Table 16

RR_no_int_no_fe<-lm(led_homs_std~product_change_export_value_share+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                      sitc_eci_2+v2csanmvch_12+v2x_rule+v2xnp_regcorr, cluster=~Country, data=xc_data)

summary(RR_no_int_no_fe)

RR_no_int_no_fe_base<-lm(led_homs_std~product_change_export_value_share+
                           v2csanmvch_12, cluster=~Country, data=xc_data)

summary(RR_no_int_no_fe_base)

#reporting
modelsummary(list(RR_no_int_no_fe_base,RR_no_int_no_fe), fmt=3, stars=T, output="latex")


#####Change in Export Value Share, Criminal Threat, and Homicides, Fixed Effects Included - main_fe
#Table 17

RR_main_model<-feols(led_homs_std~product_change_export_value_share+pop_std+gdppc_std+Ag_Emp_Percent+conflict_dummy+
                       sitc_eci_2+v2csanmvch_12+v2csanmvch_12*product_change_export_value_share+v2x_rule+v2xnp_regcorr| year+Country, cluster=~Country, data=xc_data)

summary(RR_main_model)

RR_main_model_base<-feols(led_homs_std~product_change_export_value_share+
                            v2csanmvch_12+v2csanmvch_12*product_change_export_value_share|year+Country, cluster=~Country, data=xc_data)

summary(RR_main_model_base)

#reporting
modelsummary(list(RR_main_model_base, RR_main_model), fmt=3, stars=T, output="latex")


#######Homicides, Criminal Threat and Shocks to Export Value Share - shock
#Table 19

#base model
test_3_vdem_base<-lm(led_homs_std~shock_ev_share_any+v2csanmvch_12*shock_ev_share_any, data=xc_data)

summary(test_3_vdem_base)


test_3_vdem<-lm(led_homs_std~shock_ev_share_any+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                  sitc_eci_2+v2csanmvch_12+v2csanmvch_12*shock_ev_share_any+v2x_rule+v2xnp_regcorr, data=g16)

summary(test_3_vdem)

stargazer(test_3_vdem_base, test_3_vdem)


##############prediction_plot_vdem_count.pdf - predicted homicides 
##Figure 15

test_3_vdem_no_standardized<-lm(homs_led~shock_ev_share_any+Population+GDPPC+Ag_Emp_Percent+conflict_dummy+
                                  sitc_eci_2+v2csanmvch_12+v2csanmvch_12*shock_ev_share_any+v2x_rule+v2xnp_regcorr, data=xc_data)

summary(test_3_vdem_no_standardized)


x_vec_n<-seq(min(na.omit(xc_data$v2csanmvch_12)), max(na.omit(xc_data$v2csanmvch_12)), length.out=5000)


new_dat_n_1<-data.frame(0, median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                        median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_n, 
                        median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))


new_dat_n_2<-data.frame(1, median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                        median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vec_n, 
                        median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr)))

colnames(new_dat_n_1)<-c("shock_ev_share_any", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", "v2csanmvch_12", "v2x_rule", "v2xnp_regcorr")
colnames(new_dat_n_2)<-c("shock_ev_share_any", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", "v2csanmvch_12", "v2x_rule", "v2xnp_regcorr")

fit_n_1.2<-predict(test_3_vdem_no_standardized, newdata=new_dat_n_1, interval="confidence")

fit_n_1.3<-predict(test_3_vdem_no_standardized, newdata=new_dat_n_2, interval="confidence")

plot_dat_n_1.2<-data.frame(fit_n_1.2, fit_n_1.3, new_dat_n_1)

pdf("/prediction_plot_vdem_count.pdf")
ggplot(data=plot_dat_n_1.2)+ 
  labs(x="Criminal Threat", y="Predicted Homicides")+
  geom_line(aes(x=v2csanmvch_12, y=fit_n_1.2[,1], color="No Shock to Export Value Share"),linetype="solid")+
  geom_line(aes(x=v2csanmvch_12, y=fit_n_1.3[,1], color="Shock, Export Value Share"), linetype="solid")+ 
  geom_ribbon(aes(x=v2csanmvch_12, ymin=fit_n_1.2[,2], ymax=fit_n_1.2[,3]), alpha=.3, fill="royalblue")+
  geom_ribbon(aes(x=v2csanmvch_12, ymin=fit_n_1.3[,2], ymax=fit_n_1.3[,3]), alpha=.3, fill="springgreen1")+
  geom_hline(yintercept=0.00, linetype="dotted")+
  theme_set(theme_bw())+
  theme(legend.title=element_text(size=20), legend.position="bottom", panel.grid.minor.x=element_blank(),  panel.grid.minor.y=element_blank())+
  scale_fill_manual(values=c("royalblue","springgreen"), name="")+
  scale_color_manual(values=c("royalblue", "springgreen"), name="")
dev.off()

########prediction_plot_vdem_box.pdf
##Figure 16

shock<-c(0,1)
criminal_presence<-c(0, mean(na.omit(xc_data$v2csanmvch_12)), quantile(na.omit(xc_data$v2csanmvch_12))[4], max(na.omit(xc_data$v2csanmvch_12)))

x_vals_plot<-expand.grid(shock, criminal_presence)

new_dat_box_1<-data.frame(data.frame(x_vals_plot[1], median(na.omit(xc_data$Population)), median(na.omit(xc_data$GDPPC)), 
                                     median(na.omit(xc_data$Ag_Emp_Percent)), 0, median(na.omit(xc_data$sitc_eci_2)), x_vals_plot[2], 
                                     median(na.omit(xc_data$v2x_rule)), median(na.omit(xc_data$v2xnp_regcorr))))

colnames(new_dat_box_1)<-c("shock_ev_share_any", "Population", "GDPPC", "Ag_Emp_Percent", "conflict_dummy", "sitc_eci_2", "v2csanmvch_12", "v2x_rule", "v2xnp_regcorr")
                    
fit_box_1<-predict(test_3_vdem_no_standardized, newdata=new_dat_box_1, interval="confidence", se=TRUE)

plot_dat_box<-data.frame(fit_box_1, new_dat_box_1)

pdf("/prediction_plot_vdem_box.pdf")
ggplot(data=plot_dat_box, aes(x=as.factor(v2csanmvch_12), y=plot_dat_box[,1], colour=as.factor(shock_ev_share_any), group=as.factor(shock_ev_share_any)))+ 
labs(x="Criminal Threat", y="Predicted Homicides")+
geom_point(aes(x=as.factor(v2csanmvch_12), y=plot_dat_box[,1]), na.rm=T, position=position_dodge(width=0.75))+
geom_errorbar(aes(x=as.factor(v2csanmvch_12), ymin=plot_dat_box[,2], ymax=plot_dat_box[,3]), position=position_dodge(width=0.75), width=0)+
geom_hline(yintercept=0.00, linetype="dotted")+
theme_set(theme_bw())+
theme(legend.title=element_text(size=20), legend.position="bottom", panel.grid.minor.x=element_blank(),  panel.grid.minor.y=element_blank())+
scale_color_manual(values=c("royalblue", "red3"), name="", labels=c("No Shock", "Shock to Export Value Share"))+
scale_x_discrete(labels=c("None", "Median Threat", "3rd Quartile Threat", "Maximum Threat"))
dev.off()
                          

######Drug Producing Countries and Homicides, Shock to Export Value Share - drugs
#Table 20

test_product_level_drugs<-lm(led_homs_std~shock_ev_share_any+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                               sitc_eci_2+drugs+drugs*shock_ev_share_any+v2x_rule+v2xnp_regcorr+as.factor(year)+as.factor(Country), data=xc_data)

summary(test_product_level_drugs)

test_product_level_drugs_3<-lm(led_homs_std~shock_ev_share_any+drugs+drugs*shock_ev_share_any+as.factor(year)+as.factor(Country), data=xc_data)

summary(test_product_level_drugs_3)

#reporting
stargazer(test_product_level_drugs, test_product_level_drugs_3, omit=c("year", "Country")) 

######Change in Product Export Price, Criminal Threat, and Homicides - xc_change_price
#Table 21

RR_model_change_price_no_fe<-lm(led_homs_std~diff_ex+pop_std+gdppc_std+Ag_Emp_Percent_std+conflict_dummy+
                                  sitc_eci_2+v2csanmvch_12+v2csanmvch_12*diff_ex+v2x_rule+v2xnp_regcorr, cluster=~Country, data=xc_data)

summary(RR_model_change_price_no_fe)

RR_model_change_price_no_fe_base<-lm(led_homs_std~diff_ex+v2csanmvch_12+v2csanmvch_12*diff_ex, cluster=~Country, data=xc_data)

summary(RR_model_change_price_no_fe_base)

#reporting
modelsummary(list(RR_model_change_price_no_fe_base, RR_model_change_price_no_fe), fmt=3, stars=T, output="latex")





















